      subroutine mshset
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M, ONLY:
      use cophys_com_M
      use blcom_com_M
      use objects_com_M, ONLY: chii_obj, theta1_obj, theta2_obj, nrg_obj, chi
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer , dimension(4) :: lioinput, lioprint
      integer :: n, is, it, nr, ns, kvacx, i1jk, l, ijkr, ipjkr, ijkb, ijpkb, &
                 np, kr, ix, iz, ij
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5, bbb, vvv
      real(double), dimension(idx + 2) :: thetao, thetan
      real(double) :: rhmn, pkmn, ws, ratio, fwos, fwosp, dr, rfkbp2, delt, &
         dtheta, dzstr, thetam, thetap, thetal, thetar, del6, phiang, size, &
         cphi, sphi, xx, zz, bzisq, gamma, tanhgm, phi0, wsr, pbar, wsb, wsc, &
         btht, p0, rhofac, rhomin, pknmin, wsrsav, bdbsav, xc, yc, zc, zone, &
         ztwo, xtwo, xone, wsmr, fx, fxone, fxtwo, wsnr, dummy, wk, phibc, &
         rvvol, harris,flvx,flvy,flvz,bfldx,bfldy,bfldz,t_shape,efldx,efldy,efldz,&
         xpert,zpert,exp_pert
      logical :: expnd, nomore
      character :: name*80
!-----------------------------------------------
!
!
      data rhmn/ 0.05/
      data pkmn/ 0.02/
!
      bxn(:ibp2*jbp2*kbp2) = 0.
      byn(:ibp2*jbp2*kbp2) = 0.
      bzn(:ibp2*jbp2*kbp2) = 0.
      phipot(:ibp2*jbp2*kbp2) = 0.
      rho(:ibp2*jbp2*kbp2) = 0.
      jx0(:ibp2*jbp2*kbp2) = 0.
      jy0(:ibp2*jbp2*kbp2) = 0.
      jz0(:ibp2*jbp2*kbp2) = 0.
      chi(:ibp2*jbp2*kbp2) = 0.
      qdnv(:ibp2*jbp2*kbp2,:nsp) = 0.
!
!     package for a three-dimensional torus
!     theta increases with i
!     phiang increases with j
!     r increases with k
!     rwall is the minor radius
!     rmaj is the major radius
!     scyllac package
      a = rvac
      j = 1
      if (jbp2 > 0) then
         rwz(:jbp2) = rwall
         j = jbp2 + 1
      endif
!
!ll   1.  geometric radial zoning:
      ws = rwall/rvac
      if (abs(rvac*fkbar/(rwall*float(kplas)) - 1.) > 1.E-06) then
         ratio = 1.2
         do it = 1, 35
            fwos = (ratio**kbar - 1.)/(ratio**kvm2 - 1.) - ws
            fwosp = (float(kbar - kvm2)*ratio**(kbar + kvm2) - fkbar*ratio**&
               kbar + float(kvm2)*ratio**kvm2)/(ratio*(ratio**kvm2 - 1.)**2)
            ratio = ratio - fwos/fwosp
            if (abs(fwos/fwosp) >= 1.E-05) cycle
            exit
         end do
         dr = rvac*(ratio - 1.)/(ratio**kvm2 - 1.)
      else
         ratio = 1.0
         dr = rwall*rfkbar
      endif
      rat = ratio
!
!     define rotation, iota
!
      rfkbp2 = 1./float(kbp2)
!
      k = 1
      if (kbp2 > 0) then
!      iota(k)=(k-2)*rfkbp2
         iota(:kbp2) = 0
         k = kbp2 + 1
      endif
!
!
!     define rotation matrix for toroidal geometry
!
      call rotate (ibar, jbar, kbar, rmaj, dz, strait, iota, istep, delt, &
         dtheta, dzstr, dphi, cdlt, sdlt, cdlhf, sdlhf, cdph, sdph, cdphhf, &
         sdphhf)
!
      if (nrg_obj > 0) then
!
!     vary spacing in theta to resolve a belt limiter
!
         do nr = 1, nrg_obj
            do i = 1, ibar + 2
!
               thetao(i) = (i - 2)*dtheta
               thetan(i) = thetao(i)
!
            end do
         end do
!
         do it = 1, 10
!
            wgrid(:ibar+2) = 1.
!
            do nr = 1, nrg_obj
               do i = 1, ibar + 2
!
                  thetam = 0.5*(thetao(i-1)+thetao(i))
                  thetap = 0.5*(thetao(i)+thetao(i+1))
!
                  thetal = max(theta1_obj(nr),thetam)
                  thetar = min(theta2_obj(nr),thetap)
!
                  if (thetar - thetal <= 0.) cycle
!
                  wgrid(i) = wgrid(i) + (thetar - thetal)*chii_obj(nr)/(thetap&
                      - thetam)
               end do
            end do
!
!     smooth the weight function
!
            do ns = 1, 10
               wgrid(1) = wgrid(ibp1)
               a11(2:ibp1) = 0.5*(wgrid(3:ibp1+1)+wgrid(:ibp1-1)) - wgrid(2:&
                  ibp1)
!
               a11(1) = a11(ibp1)
               a11(ibp2) = a11(2)
!
               i = 1
               if (ibp2 > 0) then
                  wgrid(:ibp2) = wgrid(:ibp2) + 0.25*a11(:ibp2)
                  i = ibp2 + 1
               endif
!
            end do
!
            call adaptgrid (3, ibar + 1, a11, a12, a13, a21, a22, a23, wgrid, &
               thetao, thetan, 0., 0., 0., 2.*pi)
!
            call gridinit_belt (ibar, jbar, kbar, iwid, jwid, kwid, delt, dphi&
               , dtheta, dx, dy, dz, rwall, rmaj, dzstr, istep, del1, del2, &
               del3, del4, del5, del6, del7, x, y, z, thetan, cartesian, xl, xr&
               , yb, yt, ze, zf, nu_len, nu_comm)
!
            i = 1
            if (ibar + 2 > 0) then
               thetao(:ibar+2) = thetan(:ibar+2)
               i = ibar + 3
            endif
!
         end do
!
      else
!
         call gridinit (ibar, jbar, kbar, iwid, jwid, kwid, delt, dphi, dtheta&
            , dx, dy, dz, rwall, rmaj, dzstr, istep, del1, del2, del3, del4, &
            del5, del6, del7, x, y, z, cartesian, xl, xr, yb, yt, ze, zf, &
            nu_len, nu_comm)
!
      endif
!
      kvacx = kvac
      if (touch) then
         kplas = kbar
         rvac = rwall
      endif
      kvac = kplas + 2
      kvp = kvac + 1
      kvm = kvac - 1
      kvm2 = kvac - 2
!
!ll   2.  set boundaries of expanded view into k=1 plane:

      if (.not.cartesian) then

         i1jk = 1
         phiang = -dphi
         size = 2.*rvac
         do j = 1, jbp2
            ijk = i1jk
            cphi = cos(phiang)
            sphi = sin(phiang)
            xx = rmaj + size
            zz = -size
            do l = 1, 4
               x(ijk) = xx*cphi
               y(ijk) = xx*sphi
               z(ijk) = zz
               go to (91,92,93,2002) l
   91          continue
               zz = size
               go to 2002
   92          continue
               xx = rmaj - size
               go to 2002
   93          continue
               zz = -size
 2002          continue
               ijk = ijk + nay
            end do
            phiang = phiang + dphi
            i1jk = i1jk + jwid
         end do

      endif

!
!ll   3.  calculate geometric coefficients:
      call geom (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibar, jbar, &
         kbar, nvtxkm, cartesian, periodicx, cdlt, sdlt, dz, x, y, z, c1x, c2x&
         , c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
         , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, vvol_half, &
         ijktmp2, a11, a12, a13)
!
!
!
!      call metric(x,y,z,tsix,tsiy,tsiz,etax,etay,etaz,nux,nuy,nuz)
      call metric (ncells, ijkcell, iwid, jwid, kwid, x, y, z, tsix, tsiy, tsiz&
         , etax, etay, etaz, nux, nuy, nuz, sie)
!
      call metricc (ncells, ijkcell, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, &
         c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, &
         c7z, c8z, vol, tsix, tsiy, tsiz, etax, etay, etaz, nux, nuy, nuz)
!
!ll   4.  calculate ic for magnetic fields:
!
!     rigid rotor screw pinch equilibrium
!     beta = total beta w/o poloidal field
!     siei is the specific internal energy on axis
!     floor on rho is rhmn*rho on axis
!     floor on p is pkmn*p0
!
      bzisq = btorus*btorus
      gamma = 0.5*log(((2. - beta) + 2.*sqrt(1. - beta))/beta)
      tanhgm = tanh(gamma)
      phi0 = rq0*btorus*rpl**2*tanhgm/((rmaj + strait*dz/2./pi)*(1. + tanhgm))
!
!     compute delta pbar
!
      wsr = -0.5*dr
      pbar = 0.
      wsr = wsr + dr
      wsb = (wsr/rpl)**2 + gamma
      wsc = tanh(wsb)
      btht = phi0/wsr*(wsc - tanhgm)/(1. - tanhgm)
      pbar = pbar - dr*btht*btht/wsr
      do while(wsc < 1. - 1.E-4)
         wsr = wsr + dr
         wsb = (wsr/rpl)**2 + gamma
         wsc = tanh(wsb)
         btht = phi0/wsr*(wsc - tanhgm)/(1. - tanhgm)
         pbar = pbar - dr*btht*btht/wsr
      end do
      wsr = wsr + 0.5*dr
      pbar = pbar - 2.*(phi0/wsr)**2
!
!     p0 is the pressure on axis
!
      p0 = 0.5*bzisq*beta - pbar
      pbar = 0.5*bzisq - pbar
      rho0 = rhoi
      if (siei /= 0.) rho0 = p0/(gm1*siei)
      rhofac = rho0/beta
      rhomin = rhmn*rho0
      pknmin = pkmn*p0
      wsrsav = 0.
      bdbsav = 0.
      inx = 2

      do k = 2, kvm
         do j = 2, jbp1
            do i = 2, ibp1
!
               ijk = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
               ipjk = ijk + iwid
               ipjpk = ijk + iwid + jwid
               ijpk = ijk + jwid
!
               ijkp = ijk + kwid
               ipjkp = ijk + iwid + kwid
               ijpkp = ijk + jwid + kwid
               ipjpkp = ijk + iwid + jwid + kwid
!
               xc = 0.125d0*(x(ipjk)+x(ipjpk)+x(ijpk)+x(ijk)+x(ipjkp)+x(ipjpkp)+x&
                  (ijpkp)+x(ijkp))
               yc = 0.125d0*(y(ipjk)+y(ipjpk)+y(ijpk)+y(ijk)+y(ipjkp)+y(ipjpkp)+y&
                  (ijpkp)+y(ijkp))
               zc = 0.125d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk)+z(ipjkp)+z(ipjpkp)+z&
                  (ijpkp)+z(ijkp))
               zone = 0.25d0*(z(ipjk)+z(ipjpk)+z(ijpk)+z(ijk))
               ztwo = 0.25d0*(z(ipjkp)+z(ipjpkp)+z(ijpkp)+z(ijkp))
               xtwo = 0.25d0*(x(ipjk)+x(ipjpk)+x(ipjkp)+x(ipjpkp))
               xone = 0.25d0*(x(ijk)+x(ijpk)+x(ijkp)+x(ijpkp))
!     calculate minor radius
!         yc=yc*toroid
!
!     wsmr is the major radius
!
               wsmr = sqrt(xc**2 + yc**2)
!
               ws = toroid*wsmr - rmaj + strait*xc
!
!     wsr is the minor radius
!
               wsr = sqrt(ws**2 + zc**2)
!
!     set the weight function to be compatible with a toroidal geometry
!
               wgrid(ijk) = wsr*wsmr
!
np=1
kr=1
call initial_state(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
!
        bxn(ijk) = bfldx
        byn(ijk) = bfldy
        bzn(ijk) = bfldz
!
!     perturbazione GEM
!
              bxn(ijk) = bxn(ijk) - pert_gem*((bx_ini*pi)/(zf - ze)*cos(2.d0*pi*&
                  (xc - 0.5d0*(xr + xl))/(xr - xl))*sin(pi*(zc - 0.5d0*(zf + ze))/(&
                  zf - ze)))
!		  /cosh((yc-.5*(yt+yb))/(yt-yb)/.1)
!     & *cos(2*pi*(yc-.5*(yt+yb))/(yt-yb)))
!
!
               bzn(ijk) = bzn(ijk) + pert_gem*(bx_ini*2.*pi)/(xr - xl)*sin(2.d0*pi*(xc - 0.5d0 &
                  *(xr + xl))/(xr - xl))*cos(pi*(zc - 0.5d0*(zf + ze))/(zf - ze))
!
!		  /cosh((yc-.5d0*(yt+yb))/(yt-yb)/.1d0)
!               *cos(2*pi*(yc-.5*(yt+yb))/(yt-yb))
!
!
!     perturbazione x-point
!
if(pert_x.gt.0.) then
xpert = xc - 0.5d0*(xr + xl)
zpert = zc - 0.5d0*(zf + ze)
exp_pert=exp(-(xpert/el_rec)**2-(zpert/el_rec)**2)

bxn(ijk) = bxn(ijk) +pert_x*bx_ini*exp_pert*( &
           -cos(pi*xpert/10d0/el_rec)*cos(pi*zpert/10d0/el_rec)*2d0*zpert/el_rec&
           -cos(pi*xpert/10d0/el_rec)*sin(pi*zpert/10d0/el_rec)*pi/10d0&
           )

bzn(ijk) = bzn(ijk) +pert_x*bx_ini*exp_pert*( &
           cos(pi*xpert/10d0/el_rec)*cos(pi*zpert/10d0/el_rec)*2d0*xpert/el_rec&
           +sin(pi*xpert/10d0/el_rec)*cos(pi*zpert/10d0/el_rec)*pi/10d0&
           )
!
               wsnr = (k - 1.5)/(kbp1 - 1.)
               phipot(ijk) = 0.5*(1. - wsnr**2)*efld0
!
end if
            end do
         end do
      end do
!
!
!	Add from a 2D file in the (x,z) plane
!
if (initial_equilibrium.eq.6) then
open(15,file=nome_file_input,form='formatted',status='unknown')

do ij = 1, ibar*kbar
!
   read(15,*)ix,iz,bbb,vvv
!   write(*,*)'read  i=',ix,'   j=',iz,'   By=',bbb
   i=ix+1
   k=iz+1
!
   do j = 2, jbp1
!
      ijk = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
!
        byn(ijk) = byn(ijk)+bbb
!
!
   end do
end do
close(15)
endif
!
!
      dummy = 0.0
!
!      call torusbc (ibp2, jbp2, kbp2, cdlt, sdlt, dummy, dz, bxn, byn, bzn, &
!         periodicx)

      call boundary_b_centered(ibp2,jbp2,kbp2, dx,dy,dz, dt, &
         bxn, byn, bzn, periodicx,1d0)

!
!
!      call bcc
!
      call vtxb (2, ibp2, 2, jbp2, 2, kbp2, iwid, jwid, kwid, vol, bxn, byn, &
         bzn, bxv, byv, bzv)
!
      call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, wate&
         (1,22))
!
      wk = 0.
      phibc = 0.
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , phipot, periodicx)
!
      call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, phipot, ex, ey, ez)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, ex, ey, ez)
!
!
!     initialize the electric field
!
      do k = 1, kbp2
         do j = 1, jbp2
            do i = 1, ibp2
!
               ijk = (k - 1)*kwid + (j - 1)*jwid + (i - 1)*iwid + 1
!
               np=1
               kr=1
               call initial_state(np,kr,harris,t_shape,flvx,flvy,flvz,bfldx,bfldy,bfldz,efldx,efldy,efldz)
!
!			   Assumes E = v x B
!
               ex0(ijk) = ex_ini * efldx
               ey0(ijk) = ey_ini * efldy
               ez0(ijk) = ez_ini * efldz
      ex(ijk)=ex0(ijk)
      ey(ijk)=ey0(ijk)
      ez(ijk)=ez0(ijk)

!
!               rvvol = 8./(vol(ijk)+vol(ijk-iwid)+vol(ijk-iwid-jwid)+vol(ijk-&
!                  jwid)+vol(ijk-kwid)+vol(ijk-iwid-kwid)+vol(ijk-iwid-jwid-kwid&
!                  )+vol(ijk-jwid-kwid))
!
!
!               ex(ijk) = ex(ijk)*rvvol
!               ey(ijk) = ey(ijk)*rvvol
!               ez(ijk) = ez(ijk)*rvvol

            end do
         end do
      end do

       exavp=ex
       eyavp=ey
       ezavp=ez
!
!ll   5.  calculate rho and csq:
      do k = 2, kbp2
         ijk = 1 + (k - 1)*kwid
         do j = 1, jbp2
            ijkr = ijk + iper
            ipjk = ijk + nay
            ipjkr = ipjk + iper
            rho(ijk) = rho(ijkr)
            vol(ijk) = vol(ijkr)
            rho(ipjkr) = rho(ipjk)
            vol(ipjkr) = vol(ipjk)
            ijk = ijk + jwid
         end do
         ijk = 1 + (k - 1)*kwid
         do i = 1, ibp2
            ijkb = ijk + jper
            ijpk = ijk + jwid
            ijpkb = ijpk + jper
            rho(ijk) = rho(ijkb)
            vol(ijk) = vol(ijkb)
            rho(ijpkb) = rho(ijpk)
            vol(ijpkb) = vol(ijpk)
            ijk = ijk + nay
         end do
      end do
!
!     *********************************
!     baal3 sets boundary conditions for cell-centered variables
!     ************************************
!      call bcc
!
      return
      end subroutine mshset
